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ABSTRACT 

We study the quasilinear evolution of the one-point probabihty density functions 
(PDFs) of the smoothed density and velocity fields in a cosmological gravitating system 
beginning with Gaussian initial fluctuations. Our analytic results are based on the 
Zerdovich approximation and laminar flow. A numerical analysis extends the results into 
the multistreaming regime using the smoothed flelds of a CDM N-body simulation. We 
flnd that the PDF of velocity, both Lagrangian and Eulerian, remains Gaussian under the 
laminar Zel'dovich approximation, and it is almost indistinguishable from Gaussian in the 
simulations. The PDF of mass density deviates from a normal distribution early in the 
quasilinear regime and it develops a shape remarkably similar to a lognormal distribution 
with one parameter, the rms density fluctuation cr. Applying these results to currently 
available data we flnd that the PDFs of the velocity and density flelds, as recovered by the 
POTENT procedure from observed velocities assuming = 1, or as deduced from a redshift 
survey of iras galaxies assuming that galaxies trace mass, are consistent with Gaussian 
initial fluctuations. 

Subject headings: cosmology — dark matter — galaxies: clustering — galaxies: 
formation 
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1. INTRODUCTION 

The "standard" model for the formation of large-scale structure is based on 
gravitational instability of small initial fluctuations in the energy density. These are 
assumed to originate from quantum fluctuations that were stretched to large comoving 
scales during the inflation phase (see Efstathiou 1990 for a review). The fluctuations are 
assumed to be a random fleld, i.e. a set of random variables, one for each point in space, 
which is fully specifled by the m-point joint probability density functions (hereafter PDFs; 
cf. Monin and Yaglom 1971). The time evolution of the PDFs for m > 2 may be sensitive 
to the nature of the dark matter (e.g. being baryonic or nonbaryonic, hot or cold; cf. 
Trimble 1987 for a review). For example, the effect of the dark matter on the two-point 
correlation function (or its Fourier transform, the power spectrum), which is a moment of 
the m = 2 PDF, is well known. But the one-point PDF is not explicitly sensitive to the 
nature of the dark matter, at least in the linear regime. This makes it a useful statistic 
for relating the present fluctuations to the initial conditions independently of the nature 
of the dark matter. 

The natural choice for the density fleld is here, as in many other physical systems, 
a Gaussian random fleld, where the one-point probability distribution of the density 
fluctuation fleld, d = dp/p, is normal with zero mean: 

P(S) = - e-^'/^'^' (1) 

(27ra2)i/2 ' 

and the joint probabilities are multivariate normal distributions. (Note that we use P to 
denote the probability density, or frequency function, and not the cumulative probability.) 
In the Gaussian case all the moments are determined by the variance a and the Fourier 
components of the density fleld have random phases. In the linear regime, the density 
fluctuations and the three components of the peculiar velocity fleld, w, are related linearly: 
5 oc —V • V, so each of the velocity components is a Gaussian random fleld too. The 
one-point PDFs of density and the three velocity components are all independent. 

However, observed correlations in the distribution of galaxies and clusters on scales 
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> 20 h-^Mpc (e.g. Maddox et al. 1990; Efstathiou et al. 1990; Bahcall 1988; Olivier 
et al. 1993), and the coherence of their velocities on large scales (e.g., Lynden-Bell et al. 
1988), have motivated a consideration of non-Gaussian initial fluctuations. This is because 
these observed correlations are in excess of the predictions of the "standard" CDM model 
when normalized to fit the distribution of galaxies on smaller scales. This model assumes 
Gaussian, adiabatic initial fluctuations with a scale-free spectrum and cold dark matter 
dominating an = 1, A = Friedmann universe. 

Theoretically, it has been shown that the general inflation picture still permits a wide 
variety of non-Gaussian fluctuations within the standard gravitational instability theory 
(e.g. Linde 1990; Kofman 1991a). Non-Gaussian perturbations also arise in other scenarios, 
e.g. where the fluctuations originate from topological defects, such as cosmic strings (see 
Bertschinger 1989 for a review) or textures (Turok 1991), or from non-gravitational cosmic 
explosions (see Ostriker 1988 for a review). The statistical nature of the initial fluctuations 
is therefore a basic distinguishing feature between major competing theories. 

Several recent observations of the large-scale fields allow determinations of the PDFs 
of density and velocity based on their spatial distributions in increasingly large volumes. 
Are these fields consistent with Gaussian initial fluctuations? Can they be used to reject 
this or other hypotheses? In order to be able to answer these questions we flrst need to 
study how the PDFs evolve in time under gravity. During linear evolution, when all Fourier 
components evolve at the same rate, the PDFs do not change form. However, nonlinear 
evolution can introduce strong non-Gaussian features. 

Weinberg and Cole (1992) have compared the effects of initial non-Gaussianity with 
nonlinear evolution from Gaussian initial conditions using a series of N-body simulations. 
They found that, while nonlinear evolution produces non-Gaussian density distributions 
which may smear out the initial conditions, some features of the initial conditions are 
preserved, enabling one to use the present PDF for distinguishing certain models of 
non-Gaussian initial conditions from Gaussian ones. 
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In this paper we study the mildly nonlinear evolution of the PDFs from Gaussian 
initial conditions. In §2 we present formal expressions for the weakly nonlinear one-point 
density and velocity PDFs in terms of the initial PDFs for general initial conditions. We 
evaluate these expressions for an initial Gaussian random field and then approximate the 
evolution using the Zel'dovich formalism in the limit of laminar flow. Then, in §3, we 
use a "standard" CDM high-resolution N-body simulation to test the approximation and 
extend the results into the multistreaming regime. Our results are applied in §4 to two 
derivations of the smoothed velocity and density fields in our cosmological neighborhood: 
the POTENT analysis of the observed radial peculiar velocities of galaxies (Dekel et al. 
1990; Bertschinger et al. 1990) and the fluctuation flelds deduced from a redshift survey 
of IRAS galaxies (Strauss et al. 1992). Our analysis and results are discussed in §5. 

2. PDF EVOLUTION IN THE LAMINAR ZEL'DOVICH APPROXIMATION 

2.1. Lagrangian vs. Eulerian PDFs 

To avoid confusion, we introduce the basic concepts and our notations methodically 
and in detail. Consider a large comoving volume V {-^ oo) in the space of comoving 
positions x, which contains a total mass M. Assume that at time t = the mass is 
distributed uniformly in space among particles of identical, inflnitesimal masses dm (— > 0), 
at initial (Lagrangian) comoving positions q. Each particle is identifled from then on by 
its Lagrangian position q. Let the Eulerian position of particle q at time t be x{q,t). If 
the mapping from g to a; is one-to-one, we call the flow laminar or single-stream. If it is 
many-to-one, i.e., if more than one q arrives at the same x at a flxed time t, we call the 
flow nonlaminar or multistream. 

Let p{x, t) be the mass density at position x at t. We may treat p as a random field in 
Eulerian space: p is drawn at random from a probability density P{p) at each position. We 
can imagine an ensemble of density fields where, at each position, P{p)dp is the probability 
that p is in the range (p, p + dp) . 

Define p(g, t) over Lagrangian space to be p[x{q, t),t], the density in the x position of 
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particle q at time t. (Note that in the multistream case, more than one Lagrangian point 
q may correspond to the same density p.) This is a random field over Lagrangian space: 
p is drawn at random from a probability density Q{p) for any randomly chosen particle. 
We again have in mind an ensemble of realizations where, for each particle, Q{p)dp is the 
probability that the particle resides in a region of density in the range (p, p + dp). The 
difference between P and Q is that the probability measure is based on volume for P and 
on mass (i.e., Lagrangian volume) for Q. 

The PDFs for the one-dimensional components of the velocity, P{v) and (5(f), are 
defined analogously. The distributions of v are assumed to be isotropic, i.e., the marginal 
distributions of each component are identical, P{vx) = P{vy) = P{vz)- 

A general relation between the Lagrangian and Eulerian PDFs of density, which will 
be useful later, is 

P{p) = P- Q{p) , (2) 
P 

where p = J pP{p)dp is the mean density. One simple way to prove this is by using the 
alternative, spatial interpretation of the distributions. Assuming ergodicity in Lagrangian 
space [that Q{p) is independent of g], one can replace the distribution over the ensemble 
of random realizations at a given q with the distribution over Lagrangian space in one 
realization, so Q{p)dp is also the fraction of mass which resides in regions where p is in the 
specified range. The total mass with that property is then MQ{p)dp. Assuming ergodicity 
in Eulerian space [that P{p) is independent of x, which, by the way, follows from the 
ergodicity in Lagrangian space when the mapping between the spaces is one-to-one and 
onto], one can replace the ensemble distribution with the spatial distribution such that 
P{p)dp can also be interpreted as the fraction of volume in which p is in the specified 
range. The total volume with that property is VP{p)dp. 

Using this interpretation involving the total mass of the particles that reside in regions 
of a given density and the corresponding volume occupied by this mass, it is clear that 
MQ{p)dp = p VP{p)dp, which implies Eq. (2) with p = M/V. 
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2.2. On Multistreaming 

Our aim in this paper is to calculate the mildly-nonlinear PDFs of density and velocity 
at time t, given the distributions at an initial time. Consider first each particle on its own. 
Imagine it to be a mass element dm initially spread uniformly inside an infinitesimal 
volume d^q. Let its volume at time t be d^Xq, where the explicit mapping from Lagrangian 
to Eulerian space, x{q,t), is provided by the dynamics. Assuming mass conservation one 
can write for each element separately g{q,t)d^Xq = pd^q, so 

(3) 

where the double vertical bars denote the Jacobian determinant. We call this density 
of an individual mass element (or "a stream"), g, the "single-stream" density. It is not 
necessarily the same as the true, "total" density p used in §2.1. This is because the mapping 
x{q) is not necessarily one-to-one. The mass elements may overlap in certain places, i.e. 
particles from different g's can cross in the same x at t. The total density p at position x 
is the sum of the contributions Qi from all the streams that arrive at x at that time. For 
a laminar fiow, g = p- 

The mapping also determines the comoving "single-stream" velocity of particle q at 
time t: '&{q, t) = dx{q, t)/dt. (Note that this coordinate velocity must be multiplied by the 
expansion scale factor to give the proper peculiar velocity.) The actual "total" velocity v 
at X is the mass-weighted average of the velocities of all the particles that cross there at 
t. The difference between •& and v is equivalent to that between the velocity of a molecule 
and the fluid velocity of a gas of molecules. 

Given the PDFs of the initial fluctuations, a specific mapping x{q, t) is used below 
(§2.3 and §2.4) to compute the Lagrangian PDFs of the single-stream quantities g and d 
at time t, which we denote Q{g) and Q{'d) (where d is one component of •&). In the case of 
laminar flow, the single-stream quantities are equal to the total ones so the single-stream 
PDFs are equal to the PDFs of the total quantities, Q{p) and Q{v). The desired Eulerian 
P{p) can then be extracted using Eq. (2) and P{v) can be extracted in a similar way. 
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dx 
dq 



The analytic results of this paper are limited to this laminar-flow approximation. It 
is, however, worthwhile to continue the analysis a bit further in the general framework 
of multistreaming in order to better understand the nature of the approximation and to 
predict its range of validity. 

Consider the general case in which multistreaming may be present. Each Eulerian 
volume element d^x may then contain several Lagrangian volume elements d^Qi, with i 
labeling the streams. The Eulerian and Lagrangian volumes are related by d^x = J{qi)d^qi, 
where J is the Jacobian determinant appearing in equation (3). We can define a total 
single-stream volume by integrating over all mass, V = J J{q)d^q. In general, V > V 
because Eulerian volumes with multiple streams are counted more than once. The Eulerian 
single-stream PDF V{q) is now defined as the probability that a point randomly selected 
from V has single-stream density in the range {g, g + dg) (and similarly for d). The 
difference between P and V is that the probability measure is based on V in the first case 
and V in the second. Using ergodicity, V{q) is also the fraction of V that has single-stream 
density in the appropriate range. Using mass conservation in analogy with Equation (2), 
one then obtains 

V{g) = ^-Q{g) , (4) 

Q 

where g = f gV{g)dg is the mean "Eulerian" single-stream density, not to be confused 
with the total p — J P{p)dp when multistreaming is important. 

An indicator for the degree of multistreaming is given by the mean number of streams 
at an Eulerian point. 

This parameter can be computed from Q{g) by 

Ns= [ ^Q{g)dg , (6) 
J Q 

where the total p is determined by the initial conditions and it changes in time only due 
to the expansion of the universe, p oc with a{t) the expansion factor. Equation (6) 



follows from the normalizatfon of as a PDF, J V{g)dg = 1, combined with relation (4) 



with g replaced by p/Ng according to the definition (5). 

Normally, A^g = 1 at early times and Ns remains close to unity as long as the flow is 
mostly laminar. If the mapping x{q, t) is continuous, i.e., it is onto the Eulerian space such 
that no empty regions are formed, A^s > 1. Eventually it grows to A^s ^ 1 in the severe 
multistreaming regime. Therefore, Ng{t) can serve as an indicator for the deviation from 
laminar flow. We will estimate Ng (t) analytically below using the Zel'dovich approximation 
to provide the mapping x{q, t). 

Given the single-stream PDFs, what are the desired total PDFs in the presence of 
multistreaming? The result for the density may be written 



where SP{p) has vanishing zeroth and first moments so that P{p) maintains the proper 
normalization and has the correct mean, p. The correction 5P is induced by caustics of the 
mapping x{q). An example of this behavior is provided by the gravitational microlensing 
problem, whose mathematics corresponds to the two dimensional Zel'dovich approximation 
extended beyond caustic formation. There the probability distribution function P{A) of 
magnification A (the analogue of density here) has the caustic-induced feature 5P{A) for 
large A (Nityananda and Ostriker 1984). The same effect is expected in three-dimensional 
dynamics, whether exact or given by the Zel'dovich approximation. 

As a rough approximation we may use equation (7) with SP = 0, which should be valid 
while Ns{t) is close to unity. However, this is certainly not an exact solution in general 
and it is not even guaranteed to be a reasonable approximation. For an exact solution one 
has to sum over a combination of single-stream probabilities under the constraint that the 
single-stream densities (or velocities) sum up (or average) to the given total density (or 
velocity). The true PDF can be schematically written as 




(7) 
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where p{N) is the probabihty that there are N streams at a point (with mean Ng) and P 
is the joint probabihty density for N streams at one point to sum up to a total density 
p, which is some function of the single-stream PDFs under the constraint. This nontrivial 
calculation is beyond the scope of this paper. 

2.3. Velocity PDF in the Zel'dovich Approximation 

In any isotropic cosmology, there is a statistical symmetry between positive and 
negative peculiar velocities (in particular (v) = 0), which must persist in the nonlinear 
regime as well. Therefore, nonlinear deviations of an initially Gaussian velocity distribution 
are subtle; they are reflected only in the fourth order irreducible moment of the one-point 
distribution — the kurtosis — or in higher even moments, which characterize features such 
as the sharpness of the peak and the extent of the tail of the distribution. (Recall that 
we are considering the PDF for one of the velocity components, e.g., Vx-) A priori, even 
such deviations could be large, but we show below that they are in fact remarkably small 
in quasilinear gravitating systems. 

Let us assume that the mapping from Lagrangian space to Eulerian space is given by 
the Zel'dovich approximation (Zel'dovich 1970), 

x{q,t) = q+D{t)7P{q) , (9) 

where D{t) is a universal function of time [D{t) oc a{t) in a spatially flat, pressureless 
universe] . The comoving peculiar velocity of particle q is then 

'&{q,t)=Xg = D^ , (10) 

and the single-stream density is given by (Eq. 3) 

qM^ ^-j-, (11) 

1 + 

^ dq 

where I is the unit tensor. For the total density p{x,t) one should sum equation (11) over 
all streams q at x. 
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It is easy to see that the Lagrangian PDF of velocities, is time-invariant under 

the Zel'dovich approximation aside from a simple scaling of in time (Eq. 10). For an 
initially Gaussian PDF we thus have for any of the three components of the velocity, at 
any time t as long as the Zel'dovich approximation is valid. 



^^^^ = MM^^'P 



^2 



alit) = D'itm') ■ (12) 



2^^(t). 

Note that the time-invariance of the form of the PDF holds regardless of the form of the 
initial PDF. 

But somewhat surprising is the fact that for Gaussian initial fluctuations and under 
the Zel'dovich approximation the Eulerian PDF of single-stream velocities is in fact equal 
to the corresponding Lagrangian PDF at all times, 

n^) = Qi^) , (13) 

so it is also time- invariant! 

To prove this, return to the interpretation of ensemble distribution, and let Q{'&, g) 
be the bivariate Lagrangian probability density for and g. Then, by Eq (4), the 
corresponding joint probability density in Eulerian space is V{'&, g) — {g/ g)Q{'&, g), so 
one can write the Eulerian PDF for velocities as 

r{^)= [ -Q{^,g)dg. (14) 
J Q 

Eq. (13) follows if ■& and g are statistically independent so that the bivariate 
distribution factors into the product of univariate distributions, 

Q(^,^) = Q(^) Q.{q) ■ (15) 

In the Zel'dovich approximation (Eqs. 10 and 11) at a fixed time the velocity of a particle 
is a function of ?/> only, while the density at a particle is a function of {d'4>/dq) only. 
Thus, if t/> and its spatial derivatives are statistically independent, then Q{'d^ g) can be 
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separated as in Eq. (15). This condition is automatically met for a Gaussian random 
field: the covariance (ipidipj / dqk) vanishes by isotropy when i/> and dijj/dq are evaluated 
at the same position. Vanishing covariance implies statistical independence for a bivariate 
normal distribution. (Equation 15 could be valid for other random fields as well but this 
issue is beyond the scope of this paper; we therefore restrict the following discussion to 
initially Gaussian fields.) 

The integral over g (in Eq. 14) can now be performed for a fixed yielding by Eqs 
(4) and the normalization of V 



proving Eq. (13). Hence, since Q('J?) is time invariant under the Zel'dovich approximation, 
V{'&) is time invariant too. 

Note that the above invariance is valid only for the one-point PDF. The velocity field 
does not remain Gaussian, in the sense that the distribution of velocities at different points 
in space is not expected to remain a multivariate normal distribution. 

We conclude that, as long as orbit crossing is negligible, P{v) = V{'d) is time invariant 
[aside from the trivial time-dependence (7^{t)] and should remain Gaussian. In the case 
of multistreaming, Q{y) and P{y) are not necessarily time invariant because, for example, 
the probability for a given number of streams at a point (Eq. 8) varies with time. Also, 
too far into the multistreaming regime the Zel'dovich approximation itself breaks down. 
We test below (§3) the validity of this result in the presence of multistreaming using the 
smoothed velocity field in an N-body simulation. 



The density PDF, contrary to the velocity PDF, is strongly affected by nonlinear 
effects. For one thing, a symmetric distribution of small density fluctuations (with (5) — 0) 
develops an asymmetry because the positive fluctuations can grow to any large value while 
the negative fluctuations are limited by deflnition to 5 > — 1 (p > 0). This eventually 




(16) 



2.4. Density PDF in the Zel'dovich Approximation 
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results in a sharp drop in Q{S) and in P{6) toward 6 = —1. Another important effect is 
due to the fact that positive (5's are typically associated with collapse, and therefore tend 
to occupy smaller volumes at later times, while negative 5's typically occur in 'voids' which 
expand in time. This tends to shift P{5) from positive to negative values. Finally, the 
formation of pancakes with high densities produces an extended tail for Q{5) and P{5) at 
large (5's. 

The key of this section is the computation of the single-stream Lagrangian PDF Q{g, t) 
of an initially Gaussian density field that has evolved under the Zel'dovich mapping. The 
somewhat elaborate calculation can be summarized as follows (cf. Kofman 1991b). Define 
Q= g/p. Based on continuity, we write the Zel'dovich density (11) as 

~Q{q:t) = ki(g,t)|-i ^ v^ = {l- DM){1 - D\2){1 - D\^) , (17) 

with \i the eigenvalues of the deformation tensor —di/ji/dqj, provided as initial conditions. 
The absolute value in the denominator allows one to continue using this expression even 
after the particle has passed through a caustic (where i^i = 0). For convenience we can 
write vi as a cubic in D{t), 

iyi = l-Dii^ + D'^ii2 - D^iJ,3 , (18) 

//l = Ai + A2 + A3, //2 = A1A2 + A1A3 + A2A3, /X3 = A1A2A3 . 



The crucial input into the desired calculation is the joint probability density for the 
eigenvalues in a Gaussian field, which has been computed by Doroshkevich (1970) to be: 



55/227 

Qa(Ai,A2,A3) = - — g- (Ai - A2)(Ai - A3)(A2 - A3) exp 



-{Sul - 7.5//2) 



(19) 



with ain equaling the variance of g at some initial time when the field is Gaussian. Equation 
(19) now allows us to determine the probability density for ^ as a function of the Aj's 
through equation (17). 
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The calculation becomes easier if we replace the eigenvalues Aj by the more convenient 
variables 

Vu U2=D^H2-D^II3, Us = D^ll3, (20) 



the first of which is simply related to the desired g (Eq. 17). Given this transformation, the 
joint PDF of the new variables can be expressed in terms of the PDF of the old variables: 



J^2, J^s) = Qa(Ai, A2, A3) 



5(l/i,I/2,J^3) 



(21) 



5(Ai, A2, A3) 

Then, the PDF of ui is given simply by double integration, 

Qui{i^i) ^ J du2 J dvzQy{vi,V2,T^z) (22) 

over the appropriate range in the {1^1,1^2, 1's) space. The desired PDF of relative density, 
Q = l^il""*^) is then given by 



(23) 



What is left is to express the Aj's in term of the i^j's and the hard part is to obtain the 
limits of integration in the double integral (22). We note when inverting the transformation 
(20) that the A^'s can be evaluated by solving the cubic polynomial equation 



A^ - D-\l -1^1 + iy2)X^ + D-'^{iy2 + J^3)A - 0-^1^3 = 



(24) 



The problem of defining the range of integration in (22) is thus reduced to finding where 
all three roots of Equation (23) are real. In practice we simply set Qu = when the roots 
are not all real. 

We finally obtain after some algebra (see Appendix A for details) 

(25a) 
(256) 



/3n(s) = s5i/2Q^^Qg 



2 . ^ 1 / 54 

- (n — \ m H — arccos — 5- — 1 
3 3 \ gs"^ 
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where a{t) = D{t)ain is the standard deviation of g = g/p ai time t according to hnear 
theory [given at tin and D{t) being the growing solution between tin and t]. The 
shape of Q{g) depends on time only via the parameter a. The numerical factor is = 
9 -5^/^/477 Ri 8.007328. / Q{g)dg = 1. In the pure laminar regime A= 1. The complicated 
expression (25) for Q{p) indeed reduces to a simple Gaussian distribution when cr <S 1 (see 
Apendix B for a proof). Otherwise, the one-dimensional integral of equation (25) has to 
be performed numerically for a given a and g. 

The desired single-stream Eulerian distribution can now be evaluated by Eqs. (4) and 

(5): 

r(i) = ^QM . (26) 

In the laminar regime equations (25) and (26) give the Zel'dovich approximations for 
Q{p/p) and P{p/p). 

We can use the derived Q{g) of the Zel'dovich approximation to estimate the mean 
number of streams at each position, Ns{a), using Eq. (6), 

Ns{a)= [-Mg)dg. (27) 
J Q 

The resultant Ng as a function of a is shown in Figure 1. A^^, much like remains flat 
at unity until cr ~ 1, and then, when orbit crossing becomes severe and the Zel'dovich 

approximation breaks down in certain places, it shoots off to large values. This growth 

is supposed to be proportional to for large a for the following reason. Replace the 
integration over Q{g)dg in Eq. (27) by 

Ns = j g-^QxiXuX2A3)dXidX2dXs . (28) 

Recall from Eq. (18) that g~^ = |1 — + — /X3-D^|. A simple symmetry argument 
guarantees that Q{Xi, A2, A3) is invariant under (Ai, A2, A3) (— Ai, — A2, —A3), so (pi) = 
(a^s) = 0- The deviation from A^^ = 1 is therefore only due to the term which grows in 
time oc D^. 
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As a bonus from this analysis we can gain, for example, some insight into the formation 
of pancakes in the Zel'dovich approximation. Simulations indicate that the process of 
pancaking is rather typical (Melott and Shandarin 1989; Nusser and Dekel 1990; Kofman 
et al. 1992). Assume that there is a universal density fall off from the nearest caustic plane, 
p{x). Such a one-to-one correspondence between x and p implies P{p) dp oc dx. Since 
we find in Eq. (25) that for large densities P{p) oc p^'^, wc get by integration x oc 
i.e. p oc a;~^/^. This is indeed consistent with the density profile of caustics in developed 
pancakes (see Shandarin and Zel'dovich 1989 for a review). Kofman (1991b), Coles and 
Jones (1991) and Coles and Frenk (1991) have also made a similar point. 

The mean number of streams at a point can tell us in turn about the rate of pancake 
formation. One can write Ng = p{l) + 3p(3) + where p{N) is the probability for N 
streams at a point (as in Eq. 8), since only an odd number of streams is possible. The 
normalization, on the other hand, guaranties that p{l) +p(3) -|- ... = 1. If Ng is not too 
big, we can ignore the occurrence of 5 streams or more and estimate the probability for 
(three-stream) pancakes to occur at any Eulerian point: p{3) {Ng — l)/2. 

3. EVOLUTION OF PDF'S IN AN N-BODY SIMULATION 

To test the nonlinear effects including the effects of multistreaming, the PDFs have 
been computed in a cosmological N-body simulation. We use a particle-mesh code 
(Bertschinger and Gelb 1991) with 256^ grid cells and 128^ particles in a periodic cubic 
box of comoving size 200 h~^Mpc. The initial conditions for the simulation are a random 
Gaussian realization of the "standard" CDM spectrum (Davis et al. 1985), assuming f2 = 1 
and h = 0.5. The spectrum is normalized such that today, based on linear growth, 
(j| = (5^) in spheres of radius 8 h~^Mpc is unity. The expansion factor at that time 
is set to a = 1. Figure 2 shows the particle distribution in the simulation at a = 1 in an 
arbitrary slice of thickness 25 h~-'^Mpc. 

Continuous density and velocity fields do not exist in an N-body simulation, where 
instead one has a point process. However, continuous fields may be defined by replacing 
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each particle with a smoothing kernel. The statistical properties of the point process are 
then reflected in those of the continuous density and velocity fields, as discussed in detail 
by Scherrer & Bertschinger (1991) and Scherrer (1992). 

The desired density and velocity fields are evaluated on a 80^ grid of comovmg spacing 
2.5 h~^Mpc. The first operations are on the grid-cell scale: a trilinear interpolation 
(cloud-in-cell assignment) of mass and momentum to a grid point from the particles in 
its neighboring cells, followed by small-scale Gaussian smoothing of the mass density and 
momentum density with a Gaussian window of radius 2.5 h~-^Mpc. The velocity at a grid 
point is then defined as the ratio of momentum to mass density there. The purpose of this 
small-scale smoothing is to obtain meaningful velocity values in grid points that reside in 
a neighborhood of empty cells so that they are not spuriously assigned zero velocity. 

The density and velocity fields are then smoothed further on a larger scale using a 
spherical Gaussian window of radius Rs- The purpose of this smoothing is to reduce the 
effects of nonlinearities by diminishing the density contrasts. A range of Rs is considered 
in order to span a range of nonlinear effects. Note that smoothing after nonlinear evolution 
does not fully remove nonlinear effects, although the longest wavelengths are expected to 
evolve according to linear theory. We will see, however, that on sufficiently large scales, 
smoothing restores approximately the linear evolution of the large-scale initial conditions. 

The smoothed velocity and density fields on the grid, at different times and with 
different smoothing lengths, are used to construct the PDFs. 

Figure 3 shows the Eulerian velocity PDF in the most nonlinear case studied with 
the simulation: a = 1 and Rg = 6 h~-'^Mpc, with ay = 277 kms~^ (and as = 0.55). The 
error bars are the standard deviation of the mean in eight octants of side 100 h~^Mpc 
each. These are only rough estimates of the true statistical uncertainties because only 
eight subvolumes were used. The PDFs at earlier times and larger smoothings are very 
similar and are therefore not shown. They are all very much Gaussian, as predicted by the 
laminar Zel'dovich approximation. The apparent excess in the positive tail beyond Scr^ is 
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at least partly due to a random deviation in the initial conditions due to the limited volume 
of the box; a similar excess shows up at all times. Thus, the N-body simulation confirms 
the Zel'dovich prediction that the velocity PDF of initially Gaussian fluctuations remains 
Gaussian in the quasilinear regime, and it extends this result into the multistreaming 
regime. 

Figure 4 shows the Eulerian density PDFs at a = 1 for three different smoothing 
lengths: Rg = 18, 12 and 6 h~^Mpc, corresponding to cr = 0.11, 0.26 and 0.55 respectively. 
The symbols are the means from eight octants in the simulation and the error bars are the 
corresponding standard deviations of the mean in the eight octants. A range of a values 
could similarily be spanned by a sequence of time steps with a fixed smoothing length. The 
dependence of the PDF on a is similar but not identical (see the discussion of moments 
and Figure 5 below). We show the smoothing sequence here because it is what one can 
obtain observationally. 

The dashed lines are based on the laminar Zel'dovich approximation [V{g) from Q{g) 
of Eq. 25 and 26] for the corresponding a values. The laminar Zel'dovich approximation 
indeed provides an excellent approximation to the simulated PDF in the intermediate, 
slightly nonlinear case with a = 0.26. At u = 0.55 the Zel'dovich approximation is still 
a very good approximation out to p/p = 3 (~ 5.5(t), but it starts to overestimate the 
positive tail beyond that. The Zel'dovich power-law tail, P oc g~^, reflects the collapse of 
the highest peaks (in Ai) into caustics — caustics which are smeared out by the smoothing 
applied to the N-body simulation. 

This example shows how smoothing and nonlinear evolution do not always commute. 
If the initial density field is smoothed first, nonlinear evolution produces caustics of high 
density. However, when an unsmoothed field is evolved and then smoothed, the smoothing 
reduces the high densities. Equation (25) is appropriate in the former case but not the 
latter. Nonlinear evolution erases some memory of the initial conditions on small scales 
but it also generates small-scale structure from the collapse of long waves (Kofman et al. 
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1992). 

The solid curves are lognormal distributions with the same u, 

'{Inp- fiiY 



where fii and ai are the mean and standard deviation of Inp. They are related to the 
corresponding moments of p, p and cr, via 

pi = \np-{l/2)al = ln(l + aV/x2) (30) 

(and recall that in fact p= 1 for p). As argued by Coles and Jones (1991) and noted earlier 
by Hamilton (1985), the lognormal distribution turns out to be an excellent approximation 
to the actual density PDF. The way it fits the simulation at cr = 0.55 over the whole range 
tested, out to ~ lOu, is striking! 

A more quantitative measure of the deviations of the PDFs from Gaussian is 
provided by their third and forth irreducible moments (Figure 5). Given a set of random 
measurement of the random variable x, define as usual the mean, p= {x), and the variance, 
cr^ = {{x — pY). Then the dimensionless skewness and kurtosis relative to the mean are 
defined by 

and ,.<(^-3. (31) 

As a reference, a Gaussian distribution has s = /c = 0. 

Based on the Zel'dovich approximation, we expect that any deviation from a normal 
shape are predominantly a function of a. We therefore plot in Figure 5 the skewness and 
kurtosis of the PDFs as a function of cr. Note though that a range of a values could be 
obtained either by analyzing the system at difi^erent times or by using a variety of smoothing 
lengths. Each panel shows three curves corresponding to three difi'erent times (a), with 
the Gaussian smoothing length (/?<, in co moving h~^Mpc) varying along each curve. The 
error bars mark the standard deviation of the mean for the moments as evaluated in eight 
octants. 
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Because of the expected symmetry between positive and negative velocities, there is 
no surprise in the fact that the velocity skewness is consistent with zero; any deviation 
must be a result of the limited volume sampled. But the fact that the velocity kurtosis 
remains constant and consistent with Gaussian is a meaningful confirmation of the laminar 
Zel'dovich approximation. One can see that the apparent small deviation from zero of both 
s and k, which is probably associated with the apparent tail in Figure 3, is indeed similar at 
all times, indicating that it must be due to the finite volume sampled rather than nonlinear 
evolution. 

While the moments of velocity remain constant over the whole range of ay tested, 
the moments of density gradually deviate from Gaussian in the corresponding range of ag- 
This deviation is significant already at relatively low ag values. 

We can also see from Figure 5 that the growth of as in the simulations is very similar 
to the prediction of linear theory: cr^ oc a to within 3% in the tested range. This is saying 
that the faster nonlinear growth in high density regions is roughly compensated for by the 
slower deepening of low-density regions (limited by 5 > —1). This allows one to use in 
the Zel'dovich approximation (25) the true, observable a instead of the linear a which we 
can't measure. 

The actual values of a and Rs that boil down to a given value of a make some 
difference, which is significant (in view of the errors) only for the density kurtosis. This 
difference depends on the power spectrum used. Since observationally a can vary only 
due to the smoothing used, we limited ourselves in Figures 3 and 4 to "observing" the 
simulation only at one time while the smoothing length is varied. Assuming a universal 
galaxy biasing factor of unity (6 = a^^ = 1), we used the time step in the simulation 
o = 1. If the true biasing factor is different from unity, then a different time step in the 
simulation should have been used to resemble the present universe (e.g. a — 0.5 for 6 = 2). 
Based on Figure 5, this would not affect much the dependence of skewness on a but the 
density kurtosis would behave somewhat differently. 
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Second-order perturbation theory predicts that the ratio of density skewness to its 
standard deviation is constant, s/cr = 34/7 (Peebles 1980), which is marked by the dotted 
hne in the top-right panel of Figure 5. We can see that this is a good approximation 
for Rg > 12 h~^Mpc, independently of how a was changed (by a or by Rs). At smaller 
smoothing lengths the N-body results tend towards s/u values in the range 4 to 3, in general 
agreement with the various models discussed by Coles and Frenk (1991). In particular, 
it has been predicted based on second order theory that this ratio would depend on the 
effective logarithmic slopes of the power spectrum at the smoothing scale, n, roughly as 
s/cr = 34/7 — (n + 3) (Juszkiewicz, Bouchet and Colombi 1992, but note they assumed 
top-hat smoothing). In our case of a CDM spectrum and Gaussian smoothing lengths in 
the range 6 — 21 h~^Mpc the effective slope is in the vicinity of n ~ — 1, which indeed 
predicts s/a ^ 3. 



4. TENTATIVE COMPARISON WITH OBSERVATIONS 

Equipped with the above results concerning the PDFs of initially Gaussian fluctuations 
for a given cr, we now make first attempts to reconstruct the one-point spatial distribution 
functions of the smoothed velocity and density fields as estimated from galaxies in a finite 
volume around us. The following comparison is tentative as it is based on limited data in a 
relatively small volume. The pilot data are provided by either the early potent analysis 
of observed velocities or the analysis of the 1.9Jy redshift survey of iras galaxies. 

The POTENT analysis used the observed radial peculiar velocities of about 1000 
galaxies in a sphere of radius ~ 60 h~^Mpc about the Local Group. The observed 
radial velocities were first smoothed into a radial velocity field on a grid, in a way that 
minimizes the effects of sparse sampling and measurement errors, potent then imposes 
the requirement of potential flow, v = — V(/), which is a natural outcome of gravitational 
instability, to reconstruct the missing two components of the velocity field (Bertschinger 
and Dekel 1989; Dekel et al. 1990; Bertschinger et al. 1990). Assuming that the galaxies are 
fair tracers of the smoothed velocity field independent of their specific type, the resultant 
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velocity field is independent of galaxy "biasing". Finally, assuming a value for Q, and 
using a quasilinear approximation for the equations governing the evolution of fluctuations 
(Nusser et al. 1991), potent yielded the mass-density fluctuation fleld that has given rise 
to the peculiar velocities. The output is provided on a cubic grid of spacing 5 h~^Mpc inside 
a spherical volume of radius 60 h~^Mpc. A Monte Carlo analysis of distance-measurement 
errors provided estimates of the uncertainties of each fleld, which enables us to use for the 
reconstruction of the PDFs only points where the uncertainty is smaller than a certain 
conservative limit, and to estimate the resultant uncertainties in the PDFs. 

Because of the limited volume analyzed, and the zero-point uncertainty in the distance 
indicators, the mean density and velocity within this volume have flnite values different 
from the universal zero means: the mean density contrast is = 0.11 and the mean 
one-dimensional velocity is Hv = 223 kms"-*^. We use the fact that if a; is a Gaussian fleld 
then the conditional probability of a; in a neighborhood where the local mean is given is 
also Gaussian with a displaced mean and somewhat reduced dispersion (depending on the 
two-point correlation function; cf. Dekel 1981, the appendix). We therefore compute and 
plot the PDFs of {S — iJ,s)/as and {v — iJ,y)/ay. 

These PDFs of potent velocities and densities, assuming f2 = 1, are shown in Figure 
6. We use only points with Monte Carlo measurement errors in the three-dimensional 
velocities and in the densities smaller than 300 kms~^ and 0.3 respectively. The resultant 
effective volume corresponds to a sphere of radius 37 h~^Mpc. The error bars are the 
standard deviation of the PDF in 30 Monte Carlo noise simulations of potent. Also 
marked are the derived flrst four moments. 

Errors due to the limited volume sampled can be estimated using the CDM N-body 
simulation, because the power-spectrum deduced from potent is not significantly different 
from the standard CDM spectrum with 6=1 (Kolatt, Seljak, Bertschinger and Dekel 
1993). In Figure 7 we show the mean PDF of eight independent spheres of radius 
50 h~-'^Mpc from the simulation and the associated standard deviations are marked by 
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the error bars. Also marked are the first four moments and their standard deviations. The 
relative volume errors for potent should be larger by a factor which is roughly the square 
root of the volume ratio, i.e. 1.57. For a conservative error estimate, the volume errors 
and the measurement errors shown in figure 7 should be added in quadrature. 

We see from figure 6 that the PDFs based on potent are consistent with Gaussian 
initial fiuctuations. However, the errors are very big, making this preliminary comparison 
only marginally interesting. Things are expected to get better though because the 
rapid progress in peculiar velocity surveys allow a potent reconstruction with smaller 
uncertainties in a larger volume. 

The IRAS analysis (Strauss et al. 1990; 1992; Yahil et al. 1990) translated the 
redshift catalog of 1.9Jy iras galaxies into a uniform galaxy- density map, whose first 
approximation was given in redshift space. The predicted peculiar velocities, and the 
corresponding corrections to the galaxy positions from redshift space to configuration 
space, were reconstructed via a self-consistent iterative scheme using linear dynamical 
theory of gravity with small-scale smoothing and quasilinear corrections (Nusser et al. 
1991), after assuming linear biasing between the density fiuctuations of galaxies and mass. 
The resultant peculiar velocity depends both on Q, and on the "biasing" parameter h 
between the density fiuctuations of galaxies and mass, but the obtained density field 
depends only weakly on the assumed value of O through the correction from redshifts to 
real positions. This analysis provided estimates for the density and velocity fields within a 
sphere of radius 80 h~-^Mpc about the local group. In the present, preliminary application 
we take these fields as provided to us by the authors of the iras analysis; we do not make 
an attempt to carefully estimate the errors involved in the sampling and in the analysis 
beyond the volume errors. 

A strong correlation is found between the iras and potent fields, both featuring as 
extended structures the Great Attractor, an adjacent large void, and the Pisces part of the 
Perseus-Pisces supercluster. The comparison yields Q,^-^/b = 1.3 ± 0.7 (Dekel et al. 1993, 
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the error bar is 95% confidence level). This allows one to adopt the simple linear biasing 
assumption as a working hypothesis and associate the PDFs derived from iras galaxies 
with the desired PDFs of the underlying dynamical mass distribution. 

The distributions of the iras densities and velocities, for two different smoothing 
lengths, are shown in Figure 8. The means are again removed and the deduced first four 
moments are marked. The errors due to the finite volume should be about twice the error 
bars shown in Figures 3 and 4 from the CDM N-body simulation, or about one half of the 
errors in Figure 7. Within these uncertainties, the iras data are consistent with Gaussian 
initial fiuctuations. The ratio s/a is consistent with being a constant as a function of 
smoothing scale, in the range s/a 1.5 — 1.8 

Note that although iras galaxies are underrepresented in rich cluster cores, so that 
they are biased against high density, this has little effect on P{p) because clusters occupy 
a very small fraction of the volume. This is an advantage of the one-point Eulerian density 
distribution over other statistics that strongly weight high-density regions. 

The more recently completed 1.2Jy iras survey (Fisher 1992) allows a more reliable 
determination of the iras PDFs in a larger volume and with more quantitative error 
analysis. The moments of the density PDF from this survey indeed seem to behave as 
expected from an initially Gaussian field subject to gravity and n ~ near the smoothing 
scale, with s/a = 34/7 — {n + 3) ~ roughly constant as a function of smoothing scale, 
in the range 1 — 2 (Bouchet, Davis and Strauss 1992). The skewness measured from the 
QDOT study of counts in their l-in-6 redshift survey of iras galaxies (Saunders et al. 
1991), given their errors, is also consistent with the above measurements (G. Efstathiou, 
private communication). 

5. DISCUSSION AND CONCLUSION 

We investigated the quasilinear effects on the one-point probability density 
functions (PDFs) of initially Gaussian fiuctuation fields. The laminar Zel'dovich 
approximation provides useful analytic expressions that are confirmed and extended into 
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the multistreaming regime using a standard CDM N-body simulation. We found that the 
velocity PDF smoothed on scales > 6 h~^Mpc is hardly affected by nonlinear evolutionary 
effects while the density PDF develops a lognormal shape. 

The observed velocity and density fields, based on potent reconstruction from radial 
velocities with $7 = 1 or on an analysis of a redshift survey of iras galaxies, have PDFs that 
are apparently consistent with Gaussian initial conditions. The data used here, however, 
are still limited in volume and they still carry large errors. Noting that random errors 
can produce a spurious Gaussian PDF, we were careful to estimate the uncertainties due 
to measurement errors and to conservatively restrict ourselves to low-error regions at the 
expense of sampling a larger volume. New data are expected to allow a significantly more 
accurate determination of the PDFs. 

These results sound encouraging for the "standard" model but can we actually use 
them to reject any non-Gaussian model of interest? Besides the fact that the current 
errors are large, we are facing several fundamental limitations. For example, there is a 
general reason for the velocity field to become Gaussian under a wide range of conditions, 
even when it came from non-Gaussian initial fiuctuations (Scherrer 1992). Recall that 
the peculiar velocity at a given point is related to the net peculiar force there (directly 
proportional to it in the linear regime and in the Zel'dovich approximation), which is the 
integral of the forces from all the mass fiuctuations around it. Assume, for example, that 
the non-Gaussianity is expressed as excessively large density peaks or wells which dominate 
the large-scale force. If the characteristic separation between these structures is not larger 
than the typical range over which the force converges, then the velocity is practically a 
sum over a few independent random fiuctuations and as such, based on the central limit 
theorem, it becomes a Gaussian variable. Thus, only certain non-Gaussian models would 
show a non-Gaussian velocity PDF, so it has to be individually evaluated for each model 
before the model can be rejected based on an observed Gaussian velocity PDF. 

For instance, the Texture model cannot be rejected based on the velocity PDF 



25 



(Gooding et al. 1992). But it still remains to be seen whether the model could be rejected 
based on the density PDF. As another example, it is clear that the non-Gaussianity scale 
in certain versions of the cosmic string model is small enough for it not to have any 
noticeable trace in the velocity PDF. A string model that could probably be rejected is 
the version where the structure is dominated by ~ 40 h~^Mpc wakes formed behind long, 
well-separated strings which accrete "hot" dark matter (e.g. Brandenberger 1991). In this 
model the velocity at a point is typically determined by the nearest wake and is therefore 
expected to be very non-Gaussian. 

The discriminatory power of the density PDF is also limited for several reasons. First, 
the deviation from Gaussianity depends on the rms fluctuation of the dynamical mass 
density, u, which is deduced from the observed a of galaxies (e.g. in the case of iras) only 
via a certain assumption concerning biasing — the relation between the galaxy density 
and the underlying mass density fluctuations which is not well determined (see Dekel and 
Rees 1987 for a review). Second, the actual deviation from Gaussianity, and in particular 
its dependence on the smoothing scale, is somewhat dependent on the shape of the power 
spectrum of fluctuations. We assumed above a "standard" CDM spectrum which seems to 
fit the data used here reasonably well, but recall that it is not the Gaussian CDM model 
that one is trying to reject here. The moral is, again, that the density PDF has to be 
evaluated individually for each specific non-Gaussian model. 

The main purpose of this paper was to provide useful tools for addressing the question 
of Gaussian versus non-Gaussian initial fiuctuations. The current preliminary comparison 
of observation with theory is encouraging for the "standard" model but it is certainly 
far from being conclusive. We hope to be able to more quantitatively rule out certain 
non-Gaussian models of interest with data that are becoming available. 
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APPENDIX A: Derivation of equation (25) 



In this Appendix we derive equation (25) for the density PDF from the double integral 
(22). The main work is to find the appropriate integration limits in (22) in terms of 
{ui, 1^2,^3), for which all three roots A(z^i, 1^2,^3) of the cubic polynomial equation (24) are 
real. 



We rewrite equation (24) in the canonical form: 



X-" + AX^ +BX + C ^ 0. 



{Al) 



First, let us find the region of the {A, B, C)-space, for which the three roots of equation 
(Al) are real. To analyze the properties of the roots of this cubic equation, we need its 
discriminant 



A = l 

4 



All three roots A are real if this determinant is negative: 



A < 0. 



(A3) 



To satisfy this condition, at least the second term on the right-hand side of equation (A2) 
has to be negative, i.e. A^ — 3S > 0. Then we can write the determinant in the form 

A = ^ C - C(-) (A, B)) C - C^+^ (A, B)) , 27C(±^ (A, B) = -2A^+9AB±2{A^-3Bf/^' 

(A4) 

The condition (A3) is satisfied if 



A^>3B and C^-\A, B) < C < C^+\A, B). 



These are the conditions required for equation (Al) to have real roots. 



(A5) 



Now we find the corresponding region in the {i'i,i'2,i^3) space where all three A are 
real. Comparing equations (24) and (Al), we have 

A ={1^1-1^2- 1)D-^ , B={v2 + V3)D-^ , C = -i^sD-^ ■ {A6) 
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Substituing equations (A6) into (A5), we get two constraints in terms of {vi, 1^2, ^3)' 

9{iyi -V2 + 2)1^3 + 9{iyi -V2- l)v2 - 2(j^i -V2-lf + 2 [{v^ -V2- \f - 8(1^2 + 1/3)] > 0, 

(^7) 

and 

- 1^2 + 2)z^3 + -^2- 1)1^2 -2{vi-V2- if - 2 [{vi -V2- \f - 3(z^2 + z^s)] < 0. 

(A8) 

In addition, the expression in square brackets must be nonnegative. We denote it 
henceforth as 

= (ui -V2- \f - 2,V2 - Sz^s, (A9). 
It is also convenient to introduce 

s = i/i - 1^2 + 2. (^10) 

The constraints become simpler in terms of vi and the new variables s and t: 

|t|3-|st2 + l(s3-27zyi)>0, -|t|3-|st2^^(s^-27zyi)<0. (All) 

The intervals of t which satisfy to both of these constraints simultaneously are 

> tl{s, vi) or tl{s, vi) <t^< tl{s, vi) , {A12) 

Here the functions tn{s, vi) are defined by 

tn{s^ ui) — s ^- + cos 6*^^ , ~ 3 a'rccos(54s~'^z/i — 1) . (^13) 

There are three roots because one may add an integer multiple of 27r/3 to 9. For 
definiteness, we label the roots according to the phase 9: Q < \9i\ < 7r/3, 7r/3 < 16*21 < 27r/3, 
27r/3 < I ^3 1 < TT, so that < < tf. In addition to the constraints imposed on t by 
equation (A12), we require that tnisjiyi) be real, implying 

s > 3^/^/^ ui > 0, (AU) 
29 



and 

s<-3|^^l|^/^ ui<0. (Alb) 

The constraints (A12)-(A15) give us the hmits in terms of t and s. They are easily 
expressed in terms of 1^2 and using equations (A9) and (AlO) to get U2 = I'l — s + 2 and 

iy3it, s, lyi) = ^{s^ -t^) -s-ui + l . {AIQ) 

Now we are ready to treat the integral (22). We have from equations (19)-(21) (the 
Jacobian in eq. 21 simply removes the three A-dependent factors from eq. 19): 



= J diy2 j dvs Q„{i'i,i'2,i^3] 



5^/^27 



du2 / di/3 exp 



3(z/i - z/2 - I? 15 



+ — =-^(i^2 + J^ + 3) 



{All) 

The first integral over 1/3 has limits given by equations (A12) and (A16) and is trivial. 
With two ranges of integration (one an unbounded interval) this integral yields three 
exponentials. The integral over V2 is performed using the variable s through equation 
(AlO). Its limits are given by equations (A14) and (A15). Their are two terms, one for 
I'l > and the other for < in equation (23). The integration of s finally reduces to 
equation (25) of the main text. 

APPENDIX B: Asymptotics of the density PDF for cr < 1 

In this appendix we show that the density PDF given as the integral (25) reduces to 
the Gaussian distribution (1) in the limit of small density dispersion cr <C 1. 

To see this, we investigate the properties of the integral (25) as cr ^ 1. The expression 
(25) can be represented as an algebraic combination of six integrals 

^ k=l "^^o 
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Here p = l/u^ » 1 is a large parameter and sq = 3g~^/^ denotes the lower limit of 
integration. The six integrals come from the six terms in equation (25) that one gets by 
multiplying out the exponentials. Four of these terms have a positive sign and two have a 
negative sign. The formulas for Fk{s) are rather tedious; some that we will use are given 
by equations (B3) and (B4) below. 

The asymptotic expansions of the integrals in equation (Bl) for p — > oo are 

Ik{p) ^ e"^^*^*-*"^ • (asymptotic series in powers of p~^) . (-^2) 

The asymptotic series have to be defined based on the analytic properties of the functions 
Ffc(s). Using equation (25b) for the functions (3n{s), after straightforward but tedious 
analysis of the functions Fk{s) in the vicinity of sq, one can show that two following terms 
in the sum (Bl) are of leading order as p — > oo: the positive term we denote I2 with 

F2{s) = {s-3f/2 + Pi{s)/2, (S3) 

and the negative term we denote with 

F3{s) = {s-3f/2 + Pi{s)/2. {B4) 

The other integrals Ii, I4, I5, and Iq are exponentially suppressed compared with I2 and 
73. Thus we have to consider the combination {I2 — h) for p — > 00. The decompositions 
of the functions ^2,3 in the vicinity of sq are: 

F.M = ^-^^ + (^) 5^ - sof' - i . (. - + .... (B5) 

F3(.) = + (^) + 5^ - .0)3/^ - i . (. - + .... (B6) 

To obtain an asymptotic series (B2), the functions Fk involved in the integrals Ik, 
have to be represented as analytic functions of the variable of integration. However, the 
series given by (B5) and (B6) show that -^2(5) and -^3(5) are not analytic functions of s. 
Therefore we must change variables of integration. For this purpose we introduce 

V = y/s — Sq. {B7) 
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Now the integrals read as 

/■OO 

Ik{p)=2 / dvve-P^'^^''\ {B8) 
Jo 

where the functions of the new variable have the following analytic series in the vicinity of 
v = 0: 

(so -3)2 /7so-3\ 2,5^ 3 1 4 



The upper sign corresponds to A; = 3 and the lower to A; = 2. For these analytic functions 
we can use the general formula of the asyptopic expansion (e.g. Olver 1974) 

Up) 2e-A(") ^^^r (!i±^) (bio) 

where T{x) is the Gamma function. The coefficients an'^^ are expressed through the 
coefficients of the series (B9). The ffist two of them, which we will use, are 

4''^) = (7so - 6)-' , (511a) 

and 

In the sum (Bl) we have the combination (I2 — Is)- Substituing in this combination 
the asymptotic series (BIO) and using (Bll), we find that the first terms with coefficients 
oo cancel. The leading remaining terms with coefficients a^'^^ give 

where we have replaced p by the original the small parameter a. 

Now we are ready to reproduce the result of the linear theory. Let us recall that 
So = Sq~^/^ fa 3 — D{t)5 for small density fluctuations 5 and a = D{t)ain- Substituing 
these formulas into (B12), we obtain 

/2 - Ig = f 9 • 5^/ V^tt) -j^ g-^'/^^*" . {BIZ) 

V / V27rcrin 

Substituing this expression into equation (Bl) we get the normal distribution for the 
density fluctuations 5. 
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FIGURE CAPTIONS 

Figure 1 : The mean number of streams at a point according to the Zel'dovich 
approximation, Ng, as functions of the hnear standard deviation of density fluctuations, 
as = a{t)ain- 

Figure 2 : The projected distribution of matter in a shce of the standard CDM N-body 
simulation used in this paper at time corresponding to linear erg = 1- The box side is 
200 h~^Mpc and the slice thickness is 25 h~^Mpc. 

Figure 3 : Velocity PDF P{v/ay) in the N-body simulation (at a = ag = 1). The 
distribution has been assembled from 80^ cubic grid points inside the box of side 
200 h-^Mpc. The error bars are the standard deviations of P in the eight octants 
of the volume. Shown is the most nonlinear case, a = 1 and Rs = 6 h-^Mpc, m 
comparison with a Gaussian distribution. 

Figure 4 : Density PDF P{p/p) in the N-body simulation (as in the previous figure, a = 
1) for three different smoothing scales: Rg = 18, 12,6 h~^Mpc, corresponding to u = 
0.11, 0.26, 0.55. The solid curves are lognormal distributions with the corresponding a 
and the dashed curves are laminar Zel'dovich distributions with the same a. 

Figure 5 : Moments of the distribution of density and velocity in the N-body simulation. 
Shown are the skewness and kurtosis as functions of the standard deviation. Each curve 
corresponds to a given time in the simulation, a = 0.5, 0.7, 1, and the variable along each 
curve is the Gaussian smoothing radius in the comoving range 6 < i?s < 21 h~^Mpc. 
The error bars are ± standard deviation in the eight octants of the volume. The dotted 
line marks the second-order prediction (Peebles 1980) ss — (34/7)a"5. 

Figure 6 : PDFs for potent velocity and density fields (il = 1) with Rg = 12 h~-^Mpc 
smoothing (Dekel et al. 1990; Bertschinger et al. 1990) in a volume (selected to 
have small errors) comparable to a sphere of radius 37 h~^Mpc. Also shown are 
the Gaussian and lognormal curves with the same a. The error bars correspond to 
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distance measurement errors as derived by Monte Carlo simulations. The relative 
errors associated with the limited volume sampled are roughly 57% larger than the 
errors shown in the next figure based on the N-body simulation. 

Figure 7 : Mean and standard deviation of the PDFs in eight disjoint spheres of radius 
50 h~^Mpc in the N-body simulation with smoothing Rg — 12 h~^Mpc. 

Figure 8 : PDFs for iras 1.9Jy density and velocity fields (Strauss et al. 1991, Yahil 
et al. 1991) in a sphere of radius 80 h'^Mpc. (a) R, = 12 h-^Mpc and (b) Rg = 
6 h~^Mpc. Also shown are the Gaussian and lognormal curves with the same a. The 
errors associated with the limited volume sampled should be roughly twice the errors 
shown in figures 3 and 4 based on the whole N-body simulation, or half the errors in 
figure 7. 
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